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ABSTRACT 

Estimating the level set of a signal from measurements is a task that 
arises in a variety of fields, including medical imaging, astronomy, 
and digital elevation mapping. Motivated by scenarios where ac- 
curate and complete measurements of the signal may not available, 
we examine here a simple procedure for estimating the level set of 
a signal from highly incomplete measurements, which may addi- 
tionally be corrupted by additive noise. The proposed procedure is 
based on box-constrained Total Variation (TV) regularization. We 
demonstrate the performance of our approach, relative to existing 
state-of-the-art techniques for level set estimation from compressive 
measurements, via several simulation examples. 

Index Terms — TV norm. Compressive sensing, FISTA 

1. INTRODUCTION 

Let X ^ W represent our signal of interest. A ^-level set of x 
is defined as the set of locations where the value of the signal x 
exceeds some specified threshold 7; i.e. 5** = *S'*(7) = {j : x{j) > 
7}, j = 1, p. Identification of level sets plays a crucial role in a 
variety of applications such as medical imaging where, for example, 
level sets can indicate presence of pathologically significant features 
such as tumors 1 1 2|. 

If X is known exactly, the level set estimation task is of course 
trivial. Here, our focus is on settings where x itself may not be di- 
rectly available; instead, we only have access to linear measurements 
of X, which may be noisy and/or incomplete. Consider, for example, 
the case where complete noisy measurements of the signal of interest 
are available, such that measurements are of the form 

y = x + n, (1) 

where y ^ M.^ and n G is an i.i.d. noise vector. In this case, a 
simple approach to level set estimation would entail coordinate- wise 
thresholding, estimating from the noisy measurements the regions 
where x exceeds the specified threshold 7. 

In more interesting settings, measurements may come in the 
form of noisy linear projections of the unknown signal, 

y ^ Ax + (2) 

where A G R'^^^ is a projection matrix and n G M'^ is addi- 
tive i.i.d. noise. Such models describe, for example, the measure- 
ments obtained in magnetic resonance imaging (MRI) applications, 
which correspond to samples of the Fourier-domain representation 
of the signal; similar models describe observations obtained in to- 
mographic imaging applications. 



Generally speaking, the condition k < p describes settings 
where the number of measurements is less than the ambient di- 
mension of the signal being acquired, a condition that may be 
imposed, for example, by sampling strategies designed to adhere 
to physical resource constraints. Problems involving recovery of 
high-dimensional signals from undersampled data comprise the 
essential focus of current research into compressive sensing (CS) 
(see, eg, 013), where more "exotic" measurement operators, such 
as (pseudo-)random projections, have been examined extensively. 
The essential goal in CS is typically to recover (or estimate) the 
unknown signal from a reduced number of measurements; here, we 
examine the problem of estimating a level set of x from compressive 
measurements. 

1.1. Prior Work 

Recent work 1 1 1 initiated the study of level set estimation from com- 
pressive measurements, and proposed a strategy for estimating the 
level set directly from the compressive measurements without per- 
forming the reconstruction step. The approach outlined in 1 1 1 com- 
prises an extension of the level set estimation technique developed 
in 1 5 1, which was designed for the setting where measurements are 
complete, but noisy, as described by the model ([TJ. 

The approach proposed and analyzed in f5l is based on a 
complexity-regularized level set estimator of the form 

S = argmin R{S) + a^{S), (3) 

where <S is a class of candidate estimates, ^(5*) is an empirical mea- 
sure of the estimator risk based on p noisy measurements of the sig- 
nal, given by 

^(^) = I - 2/^) [h^es) - Iws)] , (4) 

i 

where I^e} = 1 if event E is true and otherwise, and ^{S) is a 
carefully designed tree-based complexity regularization term which 
penalizes improbable level sets based proportional to their depth in 
a recursive dyadic representation of the signal domain. Here, a > 
is a parameter that controls the relative influence of each term in the 
overall cost function. The empirical risk function R{S) is devised as 
a surrogate for the true excess risk of a candidate level set S which 
is defined as 

v ^ — ^ 

^ 2eA(s*,s) 



where A (5'*, 5') = {i G (5'* n5')U(S'* PIS')} denotes the symmetric 
difference, and S is the complement of S. The essential idea behind 
this formulation is that estimates having small excess risk are likely 
close to the true level set, and while ^ cannot be evaluated directly 
from data (since 5** is unknown), the empirical surrogate © can. 

Now, the approach in 1 1 1 is concerned with estimating level sets 
from compressive measurements, and to that end, proceeds in a man- 
ner similar in spirit to the approach of JS), with an additional "pre- 
processing" step. Given the measured data y, the authors of 1 1 1 sug- 
gest forming the "proxy" observations 

z = A^y = X + {A^A - I)x + A^n, (6) 

" V ' 

n 

which take the form of a signal plus signal-dependent-noise h. The 
procedure of |5| is then applied to the proxy observations rather 
than y. The analysis in |1| comprises a careful treatment of the 
signal-dependent noise term n, under the framework proposed in 1 5 1. 
We refer the reader to [15] for further details and analysis of these 
existing level set estimation approaches. 

1.2. Contributions 

In this paper we examine an alternative approach, and demonstrate 
that level sets can be estimated from compressive measurements 
quickly and accurately using estimation-based techniques. Our 
method entails solving an optimization problem with total variation 
(TV) regularization, subject to additional constraints on the solution 
set. Our main contribution here is to demonstrate the effectiveness 
of this estimation-based approach to compressive level set estima- 
tion relative to a simple thresholding-based approach, as well as the 
state-of-the-art approach in 1 1 1. 

The remainder of the paper is organized as follows. In Sec- 
tion 121 we discuss the optimization problem to be solved for level 
set estimation, and briefly describe an algorithm for solving the TV 
regularization problem that is based on the Fast Iterative Shrinkage 
and Thresholding Algorithm (FISTA) We evaluate the per- 

formance of our approach via simulation, and these results are pre- 
sented in Section [3] Finally, conclusions and directions for future 
work are discussed in Section |4] 

2. TV REGULARIZATION FOR LEVEL SET ESTIMATION 

In this section, we describe an estimation-based algorithm for level 
set estimation using total variation (TV) regularization. TV regu- 
larization (proposed initially in |8 |) is now a standard approach in 
many image denoising and deblurring problems, and fast algorithms 
have been developed for solving TV minimization problems (see, 
eg., 1 7 1). Suppose we collect noisy projection measurements of (a 
vectorized version of) the image of interest, according to the model 
(|2]). As above we denote by x G the signal being acquired. Now, 
suppose p = mn for some integers m, n > 0, and let us denote by 
X the reshaped m x n image whose vectorized representation is x. 
The discrete penalized version of the TV minimization problem we 
consider here consists of solving a convex minimization problem of 
the form, 

X = argmin Haz - y\\l + a\\Z\\Tv , (J) 

Z Zi 

where A G R^^^ is the linear measurement operator, ^ G is the 
vectorized representation of Z G M'^^'', and 2/ G M^. Here, \\-\\tv 



Algorithm 1 FISTA: Level Set Estimation 

Input: L = Amax(A^A); a > 0; l,u <u) 
Initialize: /> = 1/L; = 1; x° = = 
while not converged do 

1 : X, = r'= - pVfir") 

2 : Xg = reshape(x^, m, n) 
3:Xfc=prox^(a||X||Tv)(X,) 

4 : Xfe = reshape(Xfc,p, 1) 

5 : = project(xfc, /, u) 

6 : t^^^ = (l + Vl + 4(t'^)^) /2 

7 : r'^+i = + ((t'^ - l)/t^+^) (x^ - x""-^) 
end while 



could be either the isotropic TV function, given by 

||^||TV,iso 

m — 1 n — 1 . 

m — 1 n — 1 

+ y ] \Xi^n — Xi-^l^n \ + y ] \^m,j — Xm,j-\-l\ (8) 

or the ii —based, anisotropic TV, which is defined by 

m — 1 n — 1 

||^||tv,£i = X - Xi+ij \ + \Xi^j - Xij^i\} 

i=i j=i 

m — l n — 1 

+ y ^ 1^1,71 — ^i+l,n| + y ^ \Xm,j — Xm,j-\-l\- 

The regularization parameter a > provides a tradeoff between fi- 
delity to measurements and complexity of the solution, as quantified 
by the TV norm. 

We pose our level set estimation approach in terms of the solu- 
tion of a box-constrained TV minimization 

1 2 

X = arg min - 2/II2 + A||Z||tv, 

z Z 

such that / < Zi < u, i — 1, ...,p, (9) 

where / and u are the upper and lower values every element of the 
solution X can take. In general, [l,u] may describe the range of 
allowable pixel amplitudes; for 7-level set estimation, we choose / 
to be just slightly less than 7, since for this task we are ultimately 
uninterested in those pixels whose values are less than 7. 

Algorithm 1 describes the procedure that we employ here, which 
is based on the Fast Iterative Shrinkage and Thresholding Algorithm 
FISTA [6^7}. Our optimization takes the general form 

min/(x)+^(x), (10) 

X 

where f{x) = ^\\Ax - y\\l and g{X) = a||X||Tv + Sc{X), and 
6c is the indicator function on C. The efficiency of our FISTA-based 
approach relies on us being able to quickly obtain the quantity Xk = 
prox^ (a 1 1 X 1 1 TV ) {Xg ) ; in general, for a continuous convex function 
g(x) and p > 0, 

prox^(p)(x) := arg mmlg{u) + ^\\u - x\\l \ . (11) 





(c) (d) 



Fig. 1. Sample grayscale images and their corresponding level sets. 
Panel (a) depicts an elevation map of St. Louis with pixel intensities 
in [0, 255], and panel (b) depicts its 7 = 70 level set. The image 
in panel (c) is a magnetic resonance angiography image, also having 
pixel intensities in [0, 255], and its 7 = 60 level set is shown in panel 
id). 



Here, we solve this by using the fast gradient based algorithm de- 
scribed in |7|. The function V/(x) denotes the gradient of the 
function / at the point x, which here is simply given by V/(x) = 
A^{Ax — y). The function project (x, l^u) is given by 

{X if I < X < 
I if X < /, (12) 
u if X > 1^. 

where / and u are as described above. The reshape (x, m, n) step 
simply takes a vector/matrix and reshapes it to dimension m x n. 
Standard termination criteria can be specified (eg., terminate when 
the difference between successive estimates is sufficiently small). 

3. EXPERIMENTAL RESULTS 

We performed level set estimation experiments with two different 
test imagefl shown in Fig. [Ha) and (c). Measurements of a par- 
ticular 128 X 128 test image X were obtained via the noisy linear 
model y — Ax + n, where x G (p = 128^) is the vectorized 
representation of the image X, and A G R'^^^. The images used 
for experiments are in 8-bit binary format, meaning that their pixels 
take integer values in [0, 255]. Here the entries of A are generated 
as i.i.d. realizations of A/'(0, 1/k) random variables, and the addi- 
tive noise components of y are i.i.d. A/'(0, cr^). We consider several 
settings, corresponding to several different choices of the number of 
measurements k and noise standard deviation a. 

We compare the approach outlined here (using the isotropic TV 
norm) with the state of the art approach in 1 1 1. As discussed above, 
for the proposed approach we choose / slightly smaller than the tar- 
get level (specifically, / = 7 — 5 here) and choose it = 255. In addi- 
tion, both our approach, as well as the method in (T), are dependent 

^The St. Louis elevation map image in Fig. [TJa) is available at 

www . usgs . gov/ feature s/ lewis andclark /Mapping . html; 
the magnetic resonance angiography image in Fig. [TJc) is available at 

en . wikipedia . org/ wiki /Magnet ic_resonance_angiography. 




2000 4000 6000 8000 10000 12000 14000 16000 
k 



Fig. 2. Excess risk of image from FiglHa) where the denotes 

our approach proposed here, -o- is the approach from |1| and ->- 
is the simple thresholding-based approach in a noise-free (a = 0) 
setting. 



on choice of an algorithmic parameter a > 0. Following the eval- 
uation methodology described in 1 1 1, we choose these parameters 
clairvoyantly, and display the results obtained using the parameters 
that result in the minimum excess risk, as defined in ([5]). 

FigO shows a plot of excess risk vs. number of projection mea- 
surements for estimates of the 7 = 70 level set of the image in 
Fig- da) obtained using the approach described here, the state of the 
art approach from 1 1 1, and direct thresholding of proxy observations 
defined in in a noise-free setting. The results show that the ap- 
proach outlined here achieves a much lower excess risk for the same 
number of noisy measurements. Figures [3] and IH provide a visual 
comparison of our approach with the approach of 1 1 1 when estimat- 
ing the 7 = 70 and 7 = 60 level sets, respectively, of two different 
images. Estimates obtained using each approach, for several differ- 
ent values of k (number of measurements) and different noise levels 
are shown. Here, as above, the estimates shown for each case corre- 
spond to the clairvoyantly-chosen algorithmic parameter a yielding 
minimum excess risk in each case. Note that the TV-based approach 
proposed here performs fairly well even when the number of mea- 
surements obtained are much smaller than the ambient dimension — 
see, in particular, Fig|4te) and (f). 

4. CONCLUSIONS 

We proposed a simple box-constrained total variation (TV) based op- 
timization approach for level set estimation from compressive mea- 
surements. A fast algorithm based on FISTA was discussed, and sim- 
ulation results demonstrate the effectiveness of this approach, rela- 
tive to existing techniques, in the compressive level set estimation 
problem. Future work in this direction will entail a more exhaustive 
simulation-based study of this technique for a variety of test images, 
as well as examination of the performance of this approach from 
other forms of undersampled data, such as subsampled Fourier data, 
which arises in tomographic and magnetic resonance imaging appli- 
cations. 
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